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ABSTRACT 

Several  numerical  methods  used  in  the  calculation  of  hydrodynamic 
shocks  were  investigated.  Particular  attention  was  given  to  the 
artificial  viscosity  approach  of  Von  Neumann  and  Richtmyer  and  its 
application  to  the  "PUFF"  numerical  scheme.  The  particle  model 
approach  of  Ludford,  Polachek  and  Seeger,  and  the  method  of  Lax  were 
also  considered. 
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I.   INTRODUCTION 

Frequently  the  numerical  calculation  of  comoressible  fluid  flow  is 
complicated  by  the  presence  of  shock  waves.  The  difficulties  arise 
from  the  fact  that  shock  waves  propagate  discontinuities  in  velocity, 
pressure,  and  other  variables  characterizing  the  fluid  flow.  Various 
aDDroaches  to  this  problem  have  been  offered,  each  having  its  own 
desirable  and  undesirable  characteristics.  These  approaches  generally 
fall  into  two  categories.  The  first  category  is  related  to  the  study 
of  the  viscosity  of  the  fluid,  whereas  the  second  category  is  dependent 
on  the  conservation  form  of  the  hydrodynamic  equations.  The  first 
category  of  aporoaches  will  receive  primary  attention. 


II.  NUMERICAL  CALCULATION  OF  HYDRODYNAMIC  SHOCKS 
USING  AN  ARTIFICIAL  VISCOSITY  FACTOR 

A.  INTRODUCTION 

The  equations  describing  perfect  compressible  fluid  flow,  in  the 
presence  of  shock  waves,  produce  solutions  with  discontinuities. 
Investigation  of  the  physical  situation  shows  that  the  true  discontinu- 
ities, however,  cannot  occur  due  to  the  viscosity,  or  inner  friction, 
of  the  fluid.  These  equations  ignore  the  viscosity  of  the  fluid  and 
do  not  accurately  represent  the  physical  system.  The  addition  of  vis- 
cosity terms  to  this  system  of  equations  shows  that  the  fluid  behavior 
inside  the  shock  region  is  nonlinear  but  continuous.  The  viscosity 
of  the  fluid  is  negligible  outside  of  the  shock  and  significant  inside 
the  shock.  The  original  intent  then,  was  to  replace  the  shock  reqion 
by  a  discontinuity  and  treat  the  shock  as  a  two-sided  boundary.  The 
size  of  the  jump  discontinuity,  or  rather  the  boundary  values,  would 
be  prescribed  by  the  Rankine-Hugoniot  equations.  However  this  approach 
has  several  drawbacks.  First,  the  presence  of  a  discontinuity  comoli- 
cates  the  use  of  a  numerical  scheme  to  solve  the  problem.  Secondly, 
the  shock  wave  is  in  motion,  and  hence,  the  boundary  is  movinq.  Thirdly, 
since  irreversible  thermodynamic  changes  of  state  take  Dlace  across  a 
shock  region,  an  increase  in  the  specific  entropy  of  the  fluid  must  be 
added  to  the  original  jump  conditions.  And,  lastly,  this  approach  does 
not  represent  the  physical  situation  in  that  there  is  no  indication  of 
the  behavior  of  the  fluid  inside  the  shock  region. 

Since  the  addition  of  the  true  viscosity  terms  severely  complicates 
the  system  of  differential  equations,  Von  Neumann  and  Richtmyer  [Ref.  7] 
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have  suggested  the  addition  of  a  "pseudoviscosi ty"  term  to  the  equa- 
tions of  non-viscous  flow,  which  will  not  complicate  the  system  as  much 
as  the  true  term  would.  This  term  must,  of  course,  conform  to  several 
restrictions. 

B.  THE  BASIC  EQUATIONS 

The  equations  describing  one-dimensional  flow  of  a  compressible 
fluid  are  as  follows: 

U(x.t)  =  §£  (2.2) 

p0  ||  -  -  f^Hi  (2.3) 

i^m  =  o  (2.4) 

po  at    ax  (2.5) 

e  -  $nr  (2.6, 

where  x  is  the  Lagrangean  coordinate,  X  =  X(x,t)  is  the  Eulerian 
coordinate  (i.e.,  X(x,t)  gives  the  position,  at  time  t,  of  a  fluid 
element  initially  at  position  x),  p  (x)  is  the  initial  density,  V  is 
the  specific  volume,  U  is  the  fluid  velocity,  p  is  the  static  fluid 
pressure,  E  is  the  internal  energy  per  unit  mass,  Y  is  the  ratio  of 
specific  heats  (i.e.,  c  /  c  ),  and  q  is  the  artificial  viscosity. 
Equations  (2.3),  (2.4),  and  (2.5)  are  the  equations  of  motion,  of 
energy,  and  of  continuity  respectively.  Equation  (2.6)  is  the  equation 
of  state  for  a  perfect  gas. 


C.  THE  ARTIFICIAL  VISCOSITY  FACTOR 

The  expression  for  q  must  satisfy  the  following  requirements: 

1.  Equations  (2.3),  (2.4),  and  (2.5)  must  possess  solutions 
without  discontinuities. 

2.  The  thickness  of  the  shock  must  be  everywhere  of  the  same 
order  as  ax  (the  length  increment)  used  in  the  numerical 
calculation, 

3.  The  effect  of  q  on  eqns.  (2.3)  and  (2.4)  must  be  negligible 
outside  the  shock  region. 

4.  As  ax->0,  the  solution  must  approach  a  state  with  a  jump 
discontinuity  prescribed  by  the  Rankine-Hugoniot  equations. 

Apparently,  these  requirements  are  not  enough  to  uniquely  define  q. 

In  any  case  the  expression  that  Von  Neumann  and  Richtmyer  developed  is 

,2 


q  = 


(p0cax)' 


9V 
3t 

3V 
3t 

(2.7) 


where  c  is  a  dimensionless  constant  near  unity.  By  the  use  of  equation 
(2.5)  q  can  be  written  as 


q  =  -  (cax)   3t[  . 
V    8X 


3U 


3X 


(2.8) 


Von  Neumann  and  Richtmyer  have  proven  that  q  satisfies  the  above 
requirements  for  a  particular  case  of  steady-state  plane  shock.  They 
conjecture,  however,  that  the  artificial  viscosity  approach  would  be 
equally  suited  to  more  complicated  multi -dimensional  flows.  The  problem 
they  consider  is  the  example  of  a  one-dimensional  shock  wave  separating 
two  regions  of  constant  state.  This  simulates  the  situation  that  occurs 
when  a  piston  is  pushed  at  a  constant  velocity  into  a  long  tube  con- 
taining a  fluid  initially  at  rest.  After  the  shock  has  traveled  a 
sufficient  distance  from  the  initiating  piston,  it  moves  at  a  constant 
speed,  s.  In  the  absence  of  an  artificial  viscosity  term,  the  specific 
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volume,  V,  at  some  time  t,  is  given  by  figure  1.     Since  we  are  con- 
sidering steady-state  solutions  only,  the  solutions  depend  only  on  a 
linear  combination  of  x  and  t  given  by 

w     =     x  -  st„  (2.9) 

Define 

M     =     p0S.  (2.10) 

Now 

U(x,t)  +  U(w)  =  U(x  -  st) 
implies  that  equation  (2.3)  becomes 

since 

M  Hi  =     c  1M        3W   3U         3U 
3W     P0   3W     P0  3t   3W      P0  3t 

3(p+q)  _  d(p+q)  3w    _  d(p+q) 
3x       dw     3x        dw 

Similarly,  equations  (2.4)  and  (2.5)  become 

b£+ w£  "   °  (2-12) 


and 


-«&  •  £  <2-13> 

Then  equations  (2.11)  and  (2.13)  qive 
and  equations  (2.12)  and  (2.14)  give 

d£  +  lifflffl  +  M2V  «   .    0  (2.15) 

dw  dw  dw 
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Von  Neumann  and  Richtmyer  then  integrated  equations  (2.13),  (2.14), 
and  (2.15)  with  respect  to  w  giving 

MV  +  U  =  C]  (2.16) 

M2V  +  p+q  =  C2  (2.17) 

E  +  (p+q)V  +  1/2  M2V2  =  C3  (2.18) 

as  solutions  of  equations  (2.13),  (2.14),  and  (2.15)  where  C, ,  Cp, 
and  C3  are  constants  of  integration.  Let  the  initial  and  final  values 
be  denoted  by: 

As  w-*»;  V+V.,  p^pi ,  E-*E.S  q+0  (2.19) 

As  w-»-»;  V+Vf,  p+pf,  E-+Ef,  q+0  (2.20) 

Since  V.  and  Vf  are  particular  values  of  V,  and  p.,  and  P-  are  particular 
values  of  p,  they  satisfy  equation  (2.17)  giving 

M2Vi  +  p.  +   0  =  C2 


which  implies  that 


M2Vf  +  pf  +  0  =  C2 


M2(VrVf)  =  (pf-p.)  .  (2.21) 


A  similar  argument  using  equation  (2.18)  yields 

Ef  +  pfVf  +  1/2  M2Vf2  =  C3 

-E.  -  p.V.  -  1/2  M2V.2  =  C; 


from  which  it  follows  that 


(Ef-E.)  =  1/2  M2(Vi2-Vf2)  +  p.Vi  -  pfVf  . 
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Then  by  equation   (2.21) 


<W     ■     V2  ^1  (V-2-v/)  +  p.V.   -pfVf 
=     1/2   (pf-Pi)(V.+Vf)  +  piVi   -  pfVf 
=     1/2   (p.V.   *  pfV.   -  pfYf  -  pfVf)    , 


or 


(Ef-E.)     =     1/2   (p.+Pf)(VrVf)  (2.22) 

Von  Neumann  and  Richtmyer  point  out  that  equations  (2.21)  and  (2.22) 
are  independent  of  q,  providing  q->0  as  w-*±°°,  and  in  fact  are  the  Hugoniot 
equations.  Requirement  (4)  is  then  satisfied  since  it  has  just  been 
shown  that  the  Rankine-Hugoniot  conditions  are  satisfied  for  the  flow 
sufficiently  far  from  the  shock  region.  Requirement  (4)  may  be 
examined  in  an  alternate  fashion.  Let  Z(x,t)  be  the  solution  to  the 
set  of  equations  describing  non-viscous  flow.  I.e.,  Z(x,t)  has  a  jump 
discontinuity,  prescribed  by  the  Rankine-Hugoniot  equations,  in  the  shock 
region.  Now  let  w(x,t,Ax)  be  the  solution  of  the  equations  of  viscous 
flow  corresponding  to  a  fixed  q,  or  more  accurately,  a  fixed  ax.  Then 
we  require 

limit  w(x,t,Ax)  =  w(x,t,0)  =  Z(x,t)     (2.23) 

AX-*0 

By  the  very  fashion  in  which  q  was  introduced  equation  (2.23)  is  satis- 
fied. Note  that  q  actually  has  the  dimensions  of  pressure  and  enters 
equation  (2.3)  and  (2.4)  in  the  forms   P^   and  (p+q)  respectively. 
Since  q  is  continuous,  |5~>0  as  q*0,  and  hence,  all  viscosity  terms 
approach  zero.  Consequently,  the  system  of  equations  describing  viscous 
flow  approaches  the  system  describing  nonviscous  flow  as  our  mesh  size 
approaches  zero. 

13 


The  question  naturally  arises  as  to  why  Von  Neumann  and  Richtmyer 
have  created  a  process  whereby  decreasing  the  mesh  size  causes  w(x,t,Ax) 
to  approach  a  solution  that  does  not  represent  the  physical  system. 
The  answer  is  a  heuristic  one  in  that  it  must  be  possible  to  make  q 
arbitrarily  small  to  accomodate  arbitrarily  thin  shock  waves.  It 
should  also  be  noted  that  q  is  physically  artificial  and  was  created 
for  numerical  convenience  only. 

To  investigate  the  shape  of  the  shock  Von  Neumann  and  Richtmyer 
consider  solutions  satisfying 

|^<0,  or  equivalents,  {£><><,         (2.24) 

Equations  (2.24)  are  normally  the  situation  characterizing  a  shock 
moving  to  the  right.  Then  equation  (2.7)  yields 

2 
qV  =  (Mcax)2  ({£[■)  '  .  (2.25) 

Now  from  equation  (2.18) 


E  +  1/2  M2V2  =  C3  -  (p+q)V 

E  -  1/2  M2V2  =  C3  -  (p+q)V  -  M2V2 
=  C3  -  V[(p+q)  -  M2V] 


And  then  equation  (2.17)  gives 


E  -  1/2  M2V2  =  C3  -  C2V  .  (2.26) 


Then  by  equation  (2.6) 


pV     =     (E)(Y-1) 

=   ^1m2v2  +  c3(Y-D  -  vc2(Y-i) 

B  [  ^  ]M2V2  +  C4V  +  C5  (2.27) 
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where 


C4  =  -C2  (y-D  (2.28) 


C5  =  C3  (y-D  (2.29) 

Then  equation  (2.17)  yields 

qV  =  C2V  -  pV  -  M2V2 

=  C2V  -  C4V  -  C5  -  [  2^1  M2V2  +  M2V2] 


qv 


=  C9V  -  [(4^-l  M2V2  -  C,V  -  CK  .      (2.30) 


Now  for  V  =  V.  and  V  =  V-,  q  =  0.  Therefore  since  the  right  side  of 
equation  (2.30)  is  quadratic  in  V  and  vanishes  for  V  =  V.  and  V  =  Vf 

qV     =     -[41]  M2V2  +  (C9-C.)V-C 


2     j     „      v  Vo2     o4;v      ^ 

£-    M2(-V+V.)(V-Vf)    .  (2.31) 


Y+l     »2 


Then  equation  (2.25)  yields 


(Mcax)2   [$]       .     lii    (VrV)(V-Vf) 

(<=Ax)2  [£]       =      41     (Vrv)(V-Vf)  (2.32) 

To  solve  equation  (2*32),  Von  Neumann  and  Richtmyer  proceed  as  follows 
Let, 

VVf  Vi'Vf  A 

A     =     V  -  -^-  ,     Ao     =     -Jg-1.  ,     B     =     I     (2.33) 

o 


Then, 


"*  m    -  ^1/2  <W>,/Z  <VV/2 
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Therefore, 


or 


giving 


Hence, 


where 


Finally, 


±li V2  rih     «\iai»  mV2 


=     PJ!.]         [(A0-A)(A+A0)] 


cax  £    =     P±I]1/2  [(A/-A2)]1/2   ,  (2.34) 


-'A-'  £    "     [^V2  ^V2   •  (2-35) 

0  Ao 


«*  £  ■  p£i1/2  n  -  b2]1/2  •  <2-36> 


dB 

dw  =  c!jlr]1/2cAx^(1.B2)i/a 


Y+l  (1-B<T/2 


w    =    v/    arcsin  B  ,  (2.37) 

o 


w0    =     C^]172  cax  (2.38) 


A    "    V    '    Ao  si"  5"    ■    ^  [Si"  ^L   (2-39) 

0  c  0 
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or 

VVf  ,,  VVf 

v     =    JUL    sin  £_    +    _L_L_    .  (2.40) 

2  wo  2 

V  is  obviously  continuous  and  hence  requirement  (1)  is  satisfied. 

Von  Neumann  and  Richtmyer  state  that  w  is  a  measure  of  shock  thick- 
ness. Thus,  if  c  is  near  unity,  w  is  0(ax)  and  requirement  (2)  is 
satisfied. 

Now  taking  the  derivative  of  V  with  respect  to  w  and  setting  it 
equal  to  zero  gives 

dV    VVf  1  rrtc  w     n  (9  .u 

cG7 —    w7cos  w-  "  °  '  (2-41) 

O      0 

Hence 

cos  Jj-  =  0  ,  (2.42) 

o 

giving  w  =  '  "a  '^    w  ,  (n  an  integer),  as  points  where  V  assumes 
its  relative  maxima  and  minima.  Similarly, 

A  =("^-L)1T  ■1"=-  =  0         (2.43) 
dw^       c        w  *     wo 


implies 


sin  *-    -  0 •  ,  (2.44) 

wo 


giving  w  =  n-irw  ,  (n  an  integer),  as  inflection  points  for  V.  Now 

for   W   ■    -   ye      W„     , 
c        0 


v     m     YVf  sin  (-£)  +  -i-I    =     Vf     .  (2.45) 
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For  w  =  j    w  , 


And  for  w  =  0 


VrVf  V.+Vf 

V  =  -!— I    sin  (  I  )   +  -L—L  =  v,  .     (2.46) 


v.+vf 

V  =  -!— ^  (2.47) 


Finally,  since  only  solutions  satisfying  -r-  >_  0  were  considered,  V  is 
non-oscillatory.  As  a  result  of  equations  (2.41)  through  (2.47) 
figure  2  represents  our  solution. 

Note  that  for  this  particular  problem  of  steady-state  plane  shock, 
q=0  outside  [~|  w  ,  j  w  ](i.e.,  the  shock  layer)  since  |x-  =  0  in  this 

region.  Normally,  outside  the  shock  region,  q  would  be  negligible 

2 

in  comparison  with  the  static  pressure  p  because  of  the  factor  (ax) 

in  equation  (2.7)  and  a  relatively  small  value  for  ~r.     However,  inside 

the  shock  layer  q  is  comparable  to  p  because  of  the  abnormally  large 
value  of  — p  encountered  in  that  region.  Hence  requirement  (3)  is 

satisfied.  Therefore,  for  this  particular  case  of  steady-state  plane 
shock,  Von  Neumann  and  Richtmyer  have  shown  that  their  expression  for 
q  conforms  to  all  the  necessary  restrictions. 

D.  STABILITY  OF  THE  DIFFERENTIAL  EQUATIONS 

Our  next  concern  will  be  the  effect  the  introduction  of  an  artificial 
viscosity  term  has  on  the  entire  system  of  differential  equations. 
Before  investigating  the  stability  it  should  be  noted  that  much  of 
the  computational  work  actually  done  will  be  omitted  for  the  sake  of 
brevity.  Instead,  reference  will  be  made  to  the  approach  and  ideas 
involved. 


18 


On  a  given  solution  U(x,t),  V(x,t),  etc.,  a  small  perturbation 
6ll,  6V,  etc.,  is  superimposed.  Then  the  system  will  be  stable  if  the 
perturbation  can  be  kept  arbitrarily  small  for  all  t  >  T  by  initially 
choosing  the  perturbation  at  time  T  to  be  sufficiently  small.  There- 
fore consider  the  following  variations: 

U  -  U  +  6U 

V  -  V  +  6V 

p  ■+  p  +  <sp 

q  +  q  +  6q 

We  then  obtain  our  equations  of  first  variation  (i.e.,  higher  order 

variations  are  considered  to  be  negligible): 

Pn3(gU)  =    a(6p+6q) 

3t  3X  (2.48) 


|^CT6P+(y-1)6q]  +  [Yp+(Y-l)q]  ^  +  V  ^||i 


+  |^-6V  =  0  (2.49) 


2 


*n  -  (CAX)   3U 
*     -     ~^2-     37 


9(5V)   _   3(6U) 
P0    3t        3X 


dU 


Sx 


6V  -  2%^ 


3U 


3X 


i(M  {2S0) 


(2.51) 


Equations  (2.48)  through  (2.51)  are  a  set  of  simultaneous,  linear 
differential  equations  in  fill,  fiV,  6p,  and  sq.  Their  coefficients  are 
composed  of  terms  depending  on  the  solution  functions  U,  V,  p,  and  q. 
Since  U,  V,  p,  and  q  are  considered  to  be  smooth,  well-behaved  functions 
of  x  and  t,  they  will  be  considered  to  be  constants  in  a  small  region. 
Equations  (2„48)  through  (2,51 )  were  combined  into  one  equation  and  a 
separation  of  variable  technique  was  employed.  Using  this  technique 
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possible  solutions  were  of  the  form: 


(2.52) 


<5U  =  6Uoeat  (cos  kx  +  1  sin  kx)  =  sUoelkx+at 

6V  =  6VQeat  (cos  kx  +  i  sin  kx)  =  6V0eikx+at 

6p  =  6pQeat  (cos  kx  +  i  sin  kx)  =  6poelkx+at 

6q  ■  6qQeat   (cos  kx  +  1  sin  kx)  =  6qoeikx+at 


where  SU  ,  SV  ,  fip  ,  fiq  ,  k,  and  a  are  constants  and  k  is  real.  Sub- 
stitution of  equations  (2.52)  into  equations  (2.48)  through  (2.51) 
yields  the  following  set  of  simultaneous  homogeneous  linear  equations 


in  5Uo,  6Vo,  6p0,  and  6qQ: 

RS     =     0 
where  R  is  the  following  4x4  matrix: 


(2.53) 


ap, 


ik 


E$r  +"V] 


R     = 


[2ik^ 


-ik 


9U 


8X 


ik  0 

9V- 


[(y-l)|f]        [TPa+(Y-l)qa-^] 


[- 


(CAXT    8U 


V 


2         9X 


ap 


(2.54) 


and  S  is  the  following  vector: 


ikx+atN 


i  kx+at 


i  kx+at 


i  kx+a1 


(2.55) 
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Now  if  equation  (2.53)  is  to  hold  for  a  nontrivial  S,  then 

DET(R)  =  0 
The  characteristic  equation  is: 


(2.56) 


,  sZ   ;     9poa  3V  r!/rAVN2 
(aP0)  Y3t  +  2—Jt   (kcAx) 


,2  9V 


3U 

3X 


-Jrtto-i'g 


3D 


3X 


+  k  a[yp+(Y-l)q] 


+  pQ2a3V  +  2pQa2(kCAx)2 


3U 
3X 


.  °L  (kcAX)2  3U 
V  lKCAX;   3X 


3JJ 

3X 


+  k2 15-  =o 

3t 


(2.57) 


Equation  (2.57)  establishes  the  relationship  between  a  and  k.  By  fixing 
k  and  examining  the  corresponding  a,  we  can  investigate  the  behavior  of 
the  perturbation.  It  should  be  noted  that  the  behavior  of  the  pertur- 
bation must  be  examined  both  in  the  shock  regions  and  in  the  normal 
regions.  In  the  shock  regions  all  terms  will  be  retained.  In  the 
normal  regions,  terms  containing  dissipative  factors  (i.e.,  terms  con- 
taining ax)  are  dropped.  Now  we  are  interested  only  in  the  a's  corres- 
ponding to  very  large  k's.  Hence  we  retain  only  the  dominant  terms, 

2  3 
in  a  and  k,  of  equation  (2.57).  The  dominant  term  in  a  is  p  a  V. 

2 
The  dominant  term  in  k  contains  k  and  is  therefore  dominated  by 

3U 


2p  a2(kcAx)2 

0 


3X 


,  which  is  the  dominant  term  in  ak.  Hence  in  the 


shock  regions,  equation  (2.57)  reduces  to: 


p  2a3V  +  2p  a2(kCAx)2 


3U 

3X 


=  o  , 


(2.58) 


giving 


a  = 


-2(kcAx)' 


>oV 


3U 
3X 


(2.59) 
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and  in  the  normal  regions 


giving 


,2       2  3,, 
koryp  +  pav 


2 

k  yp 


=  0, 


(2.60) 


(2.61) 


Thus,  in  the  normal  regions  our  original  system  of  differential 
equations  is  stable,  and  in  the  shock  regions,  the  system  is  asymptotically 
stable.  Von  Neumann  and  Richtmyer  also  point  out  the  terms  in  the 
equations  of  variation  that  lead  to  the  dominant  terms  fn  equation 
(2.57).  In  the  shock  regions: 


3(6U) 

3t 


32(6U) 
1   2 

3X 


where 


a  = 


2(cax)' 
Vp„ 


3U 


3X 


and  in  the  normal  regions 


32(6U)  ~  -  2  32(6U) 
3t2      °    3X2 


(2.62) 


(2.63) 


(2.64) 


where 


2  =  IP_ 


2\/ 
p0  V 


(2.65) 


E.  FINITE  DIFFERENCE  SCHEME 

To  solve  the  system  of  differential  equations  several  finite 
difference  schemes  could  be  used.  The  one  Von  Neumann  and  Richtmyer 
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offer  is  ingeniously  simple.  The  central  differences  used  are  skill- 
fully staggered,  taking  advantage  of  the  artificial  nature  of  q. 

Let  our  rectangular  mesh  with  increments  ax  and  At,  and  Integers 

i  =  0,  1,  2,  . ,.,  L;  n  =  0,  1,  2,  ...  be  contrived  in  the  following 


fashion: 


Ui+l/2  =  uC(i+V2)Ax,  nAt] 


Vi+l/2  =  V[(i-H/2)ax,  nAt]  ,  etc.        (2.67) 

After  rewriting  equation  (2.4)  the  finite  difference  scheme  is  as 

follows:  n  1/2 

Un+1/2  -  Un_1/2  pi-H/2  +  Vl/2  "  Pj-1/2 

p      _, t = 

At  AX 

n-1/2 
qi-l/2  (2.68) 

AX 


vn+l    yn       nn+l/2   ,,n+l/2 

Vl/2  "  M+1/2  ui+l      "  uf 

p    =    

0  At  AX 


(2.69) 


nn+l/2               2(cAx)2(U?:]/2-Un+1/2)|u^/2<1/2 
qi+l/2     "     '  — — 

(AX)2  {+wn  ♦  V^]/2)   (2.70) 


r.   Phi/2  +  Pi+1/2  +   (Y-l)q^]^]   (Vitl/2  [  V1+1/2> 


,.,n+1  un        wnn+l  n        % 

+     (Vi+1/2  +  Vl/2)(Pj+l/2  -  Pi-H/2}       m     0     {2n) 

2At 

Since  central  differences  were  used  the  discretization  error  will 
be  0(ax)2  and  0(At)2. 
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F.  STABILITY  OF  THE  FINITE  DIFFERENCE  SCHEME 

Having  shown  the  original  system  of  differential  equations  is 
stable  and  having  chosen  a  suitable  finite  difference  scheme  one  must 
now  show  that  the  difference  equations  are  in  fact  stable.  A  finite 
difference  scheme  is  said  to  be  unstable  if  the  rounding  error,  intro- 
duced in  approximating  the  numerical  solution,  grows  exponentially  with 
each  iteration,  making  nonsense  of  our  numerical  data.  This  concept  is 
analogous  to  the  one  for  differential  equations  in  that  the  rounding 
error  corresponds  to  a  small  perturbation  in  the  numerical  solution. 

Von  Neumann  and  Richtmyer  have  shown  equations  (2.68)  through 
(2.71)  to  be  conditionally  stable  (i.e.,  stable  only  for  certain  com- 
binations of  ax  and  At).  Hence  certain  restrictions  will  be  placed 
on  the  choice  of  ax  and  At.  It  should  be  further  noted  that  these 
restrictions  will  not  be  the  same  in  shock  regions  as  in  normal  re- 
gions .  As  before  much  of  the  computational  work  actually  done  will  be 
omitted  for  brevity  and  clarity. 

Outside  shock  regions  the  stability  criteria  of  Von  Neumann  and 
Richtmyer  is 

|i(!V-)    <1.  (2.72) 

po  V 

Equation(2.72)  is  the  usual  stability  criteria  encountered  when  hydro- 
dynamic  equations  of  the  form  of  equation  (2.64)  are  approximated  by 
central  differences. 

A  similar  analysis  in  the  shock  regions  yields 

2aAt 


(AX) 


>  £  1.  (2.73) 
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From  equations  (2.38),  (2.39),  and  (2.62) 


,  .  Scax  l^p-l   [I|l]1/2  cos  J-  .      (2.74) 


Equation  (2.73)  then  becomes 


where 


Now  let 


J+l    .   c.n  w_ 
At     ,2'/2     M  Wo  (2.75) 

Ax1y+1  4sccos^ 

wo 


0     ■    w-  (2.76) 

vf 


J+1     +    ei-n    W 

«S->  -  EI       ^ 

0  w 

COS   ~ 

wo 
=    CT^iT+^S-    -  (2'77) 

O  0 


for  -tt/2  w     <  w  <  tt/2  w  . 
Then, 


■t-y  sec  —  tan  —  +  sec  — 

J-  I     w^     w„       w„ 

0       0  0 

w_ 
J+l  sinwo    +    1 


J-T  „  2  w     „  2  w 
cos  —    cos  — - 

0  0 


-  sec2^(£j)sin^+l]        (2.78) 

0  0 


Now  setting  F'  =  0  implies 
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sin  g_  ,.  .  (£})    ,  (2.79) 

0 

since  sec  e  is  never  zero.  Hence,  we  have  obtained  as  a  critical  value 

J+l  J^ 

F     =     J-l  J+l 


[i-^4]  ^ 
(j+D2 


2J1/2 

=      ^rr  (2-80> 

Noting  that  as  j-j-  ■>  5"    from  tne  left  Ffi~)  "*"  °°-     Since  only  one  critical 
o  o 

point  has  been  found,  equation  (2.80)  must  be  the  minimum  value  for 

F(jJ-)  in  the  shock  region  (-tt/2  w    <  w  <  tt/2  w  ).     Consequently, 
o 


equation  (2.75)  may  be  rewritten  as 


Hi    ^V2    Iw]  (2-81) 


From  equation  (2.21) 


(Pf-P{)  I/O 

M  ■  [tv^t]       •  (2-82) 


Equations  (2.22)  and  (2.6)  give 


p-V 


so  that 


fVf     _     PiVi     =     1/2   (p1+Pf)(VrVf) 

=     1/2  pf(V.-Vf)  +  1/2  pf(V.-Vf), 

(2.83) 

PfVf  i  P.-Vf 

-1^  -   1/2  Pf(V.-Vf)     =     [1/2  +  lrr]p.Vi   -  -J-I  ,   (2.84) 
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yielding 


2pfV 


+1, 


£f  -Pf  (vrvf)  =  t(SJvt^  Pl 


^f 


Then  equations  (2.10)  and  (2.21)  give,  after  eliminating  p., 


S  =  —  1 


Vf  -  Pf  (y,  -vf ) 
hr-j  J 


r+1 


L(^)vrvf 


(vrvf) 


1/2 


(2.85) 


or 


S  = 


1  (y+1)-J(y-1)  -i 
7  (v+1)J-(y-D 


J-l 


V2 


Finally 


1/2 


p0  lJ(y+D-(y-D- 


LY  vf  J 


1/2 


1  r  .  .2  .  ..  i 

r0    lj(y+d-(y-d  J 


1/2 


'of 


(2.86) 


where  S  -  is  the  speed  of  sound  behind  the  shock,  relative  to  x. 
Now  equations  (2a86)  and  (2.81)  give 


At      1 


ax  -  2S  fc 


(y-D 


[  J"TFTT]J 


S/2 


(J-l)1 


(2.87) 


Since  Vf,  and  hence  J,  is  generally  unknown  until  the  problem  has  been 
solved,  the  term  inside  the  radical  is  replaced  by  its  minimum  with 
respect  to  J: 
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Let 


Then 


[  (j.ri)  J]  1/2 

Q(J)  =  3±! .         (2.88) 

(J-D 


-1/2  1/2 

1/2(J-1)[(J-^|)J]-[2J^]-[(J-^|)J] 

(J-D2 
1/2(J-1)(2J-l^-)-[J-I^]J    t    (2<8g) 
?      l     V2 

(j-D2  ccj-sJ-)  j: 


Setting  Q'  =  0  gives 


so  that, 


But 


j£T/2(^})-  1]  +  V2(^|)=  0         (2.90) 


J  =   ri  .  (2.91) 


—    —   V-  1 


r+1 


since  J  =  ^— r  corresponds  to  an  infinitely  strong  shock  and  J  <  1 

Y-l 

implies  that  Q  is  imaginary.  Since  ■L^  <   1,  there  are  no  critical 
points  in  the  interval  [1,  ^4-].  Hence  J  =  ^g  is  the  P°int  where 
Q  assumes  its  minimum  since  Q  can  be  made  arbitrarily  large  as  J+1 
from  the  right.  Hence 


<u 


K  $(£}>- 1]1/2 


mm 


Y+l 


1/2  .  (2.93) 

Y 
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Finally  substitution  of  (2.93)  into  (2.87)  yields 

At  V2 

fli-*—      •  (2.94) 

"      2sofc 

Equation  (2.94)  is  then  a  sufficient  condition  for  stability  of  the 

difference  equations  inside  the  shock  region.  It  should  be  noted  that 

Von  Neumann  and  Richtmyer  have  ignored  boundary  conditions  in  their 

stability  analysis.  Should  derivative  boundary  conditions  enter  the 

problem  an  appropriate  stability  analysis  will  have  to  include  a 

discussion  of  the  difference  equations  used  to  approximate  the  boundary 

conditions. 
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III.   INTEGRATION  OF  THE  EQUATIONS  OF  TRUE  VISCOUS  FLOW 

An  alternate  method  to  the  numerical  calculation  of  hydro dynamic 
shocks  is  offered  by  Ludford,  Polachek  and  Seeger  [Ref.  5].  Their 
approach  is  similar  to  that  of  Von  Neumann  and  Richtmyer  in  that  their 
analysis  examines  the  viscosity  of  the  fluid.  However,  by  replacing 
the  fluid  continuum  by  a  particle  model,  they  show  that  the  true  equa- 
tions of  viscous  flow  can  be  numerically  integrated. 

A.  THE  BASIC  EQUATIONS 

The  equations  of  one  dimensional  flow  of  a  perfect  viscous  compres- 
sible fluid  may  be  written 

DU      8p_  da  /-  ,n 

p  Dt  "  --^  +  ^x"  '  (3J) 

Dp  _     3U  ,,  0, 

Dt  "  "p8x"  (3'2) 

aw  -  PT^  (3.3) 

p  =  PRT  J  R  =  Cv(Y-l),  (3.4) 

where 

_  4  3U  /,  cx 

a  "  3^  33T  '  (3'5) 

S  =  Cv  log(p/py).  (3.6) 

T,  a,  S,  and  u  are  the  temperature,  viscous  stress,  entropy  per  unit 
mass,  and  coefficient  of  viscosity  respectively.  D  represents  the  rate 
of  change  of  the  quantity  written  after  it,  if  we  move  with  the  gas 
particle  (Lagrangean  viewpoint).  All  other  quantities  were  defined 
previously.  Now  p  may  be  eliminated  from  eqn.  (3.3)  by  use  of  eqns. 
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(3.4),  (3.6)  and  (3.2)  in  the  following  manner 


9U     T  DS    _p  /Tx  DS 
0   &T  =  pT  Dt  =  =f  (T)  Dt  ' 


giving 


Since 


we  have 


Then 


finally  giving 


3U  _    p    DS  /0  7N 

°^  ~  y^ry  fit  ■  (3-7) 


ft     -     Cv^(log(p/p^)}  .  Cv<  !{£-*{£>    (3.8) 


0&L    =  _£_  {  lDp__  XDp_  ,  Q) 


aW  -  VFU  W   ~    (FT)  p  97        (3J0) 


8|  -  f  [c(Y-l)-YP]  (3.1D 


Note  that  eqn.  (3.11)  is  independent  of  p. 

Now  Ludford  et  al .  approximate  the  motion  of  the  fluid  by  the  motion 
of  a  system  of  particles.  They  consider  flow  inside  a  tube  of  unit 
cross-section.  The  tube  is  considered  to  contain  (N+2)  particles,  each 
of  mass  m,  and  each  moving  along  a  fixed  straight  line.  Each  particle 
then  represents  that  mass  m  of  gas  which  initially  had  the  particle  at 
its  center  (see  figure  3).  Hence  the  forces  encountered  within  the 
original  gas  may  be  approximated  by  the  interaction  of  forces  between 
the  particles. 

■f"h 

If  x  is  the  x  coordinate  of  the  n   particle,  and  a  dot  represents 
differentiation  with  respect  to  time,  the  viscous  interaction  between 
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the  n       and  (n+1)       particles  is  given  by 


r-Vl   "  xn 


an+%    =    Iy  Ex„,1   -  x„  <3-12> 


^n+1       An 
as  a  result  of  eqn.    (3.5).     Now  from  eqn.    (3.1) 


-(Pn+>  -  Pn_jJ         (c  +p  -  o     , ) 

pxn  x  ~    -  x     ,)     +  (x  ~    -  x     ,)     '     (3J3) 

v  n+V       n-v         K  n+k        n-V 


giving 

mx     =     "(Pn^-  K-p)  +  (an^"  <W    •  ^3-14) 

since 

m    =    p(V%-V^(D-  (3-15) 

The  pressure  interaction  between  the  n      and  (n+1)       particles  is  given 


by 


x  ,  -  x 

P^  =  xn+1  -  x"  [{^K+k   -  YPn+%3-     (3-16) 
n+ 1    n 


If  the  constants  u  and  y   are  those  of  the  original  gas  equation  then 
eqn.  (3.2)  is  automatically  satisfied.  Hence  eqns.  (3.12)  through 
(3.16)  are  sufficient  to  characterize  the  particle  model. 

Before  choosing  a  finite  difference  scheme,  Ludford,  Polachek,  and 
Seeger  nondimensionalize  eqns.  (3.12),  (3.13),  (3.14),  and  (3.16)  by 

the  following  transformations: 

1     2 
P  =  P.P  =  L  P«a  ^P 
K    Ko     y  °  ° 


(3.17) 


0       = 

n       * 

X      = 

U 

t    = 

* 

£T 

V 

2 

YPo 

Pn 
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qivinq  in  olace  of  the  original   equations  of  narticle  motion, 


t  _     fcl     -  Xn 


>'<*  -   j^j-   Kt-1)*^  T^l  • 


where  primes  denote  differentiation  with  respect  to  T,  and 


A=  (j)  r^-  Cy(N-h1)]1/2  (3.19) 


B.      FINITE  DIFFERENCE  EQUATIONS 

The  finite  difference  scheme  chosen  to  represent  eqns.    (3.18)   is 
as  follows: 

Ym+1       Ym  Ym+1  Ym 

vm+^     =       2  r(  n+1   "  Vl)   -  (   n       '  Xn)-, 

n+k          AT  L/Ym+^       Ym     x  /Ym+1  vitkj 

Un+1  +  VV  "  Un      +  V 


C*  -  y&$  +  «#■  <3-20) 


X^-ZX^X-1      -      (AT)2:-(P-       -   P-     ) 


m        „  m 


+  A(Z"'      -£  '"  ,n 


"CI-Ck-  ^S3|K(r-i)/P*S-j(C.*0 


33 


m  +■  V\ 

In  the  system,  the  quantity  P  v     is  the  pressure  between  the  n 
and  (n+l)th  particles  at  time  T  =  mAT.  It  should  be  noted  that  the 
above  system  cannot  be  solved  explicitly  since  the  superscripts  of  I 
are  staggered.  Hence  eqns.  (3.20)  must  be  solved  by  an  iterative  pro- 
cedure at  e^ery     AT/2  time  increment. 

The  stability  analysis  of  the  finite  difference  equations  is  quite 
similar  to  that  of  Von  Neumann  and  Richtmyer.  Equations  (3.20)  are 
numerically  stable  if,  in  the  normal  regions, 

(Xi  -  X  ,)  1/2 

AT  <  [  "  yP  "  ]  ]    •  (3.21) 

The  system  is  unconditionally  stable  in  shock  regions. 

The  aooroach  taken  by  Ludford,  Polachek,  and  Seeger  has  some 
advantaqes  over  that  of  Von  Neumann  and  Richtmyer.  First,  in  the 
particle  model  approach,  the  stability  requirements  are  somewhat  less 
stringent.  Secondly,  this  approach  gives  a  fairly  accurate  representation 
of  the  behavior  of  the  fluid  in  the  shock  regions,  whereas  the  artificial 
viscosity  approach  does  not  represent  the  physical  system  in  the  shock 
layers.  However  the  approach  by  Ludford  et  al .  does  have  a  major  dis- 
advantage. If  the  fluid  under  investigations  has  a  very  low  viscosity 
an  inconveniently  fine  mesh  will  have  to  be  used  in  order  that  the  finite 
difference  approximations  are  accurate  beyond  the  shock  front.  Otherwise 
AX  will  be  larger  than  the  thickness  of  the  shock  wave.  It  should  be 
recalled  that  in  the  artificial  viscosity  approach  q  automatically 
adjusted  to  the  shock  thickness. 
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IV.  THE  METHOD  OF  LAX 

It  was  stated  earlier  that  approaches  to  the  numerical  calculation 
of  comoressible  fluid  flow  generally  fall  into  two  categories.  The 
first  category  consisted  of  methods  which  examined  the  viscosity  of  the 
fluid.  The  motivation  for  this  approach  was  that  when  a  nonlinear  term, 
multiplied  by  a  small  coefficient,  is  introduced  into  a  differential 
equation,  it  may  produce  large  changes  in  the  behavior  of  the  solution. 
The  second  category  consists  of  methods  which  examine  the  basic  equations 
of  nonviscous  flow,  while  allowing  discontinuous  solutions.  An  ingenious 
method  developed  by  Peter  Lax  [Ref.  4]  belongs  to  the  second  category. 

It  is  evident  that  the  approaches  discussed  so  far  are  heuristic  in 
nature  and  lack  theoretical  preliminaries.  Lax  succeeded  in  makinq 
several  theoretical  observations  which  could  tie  together  the  various 
methods  discussed  in  this  paper.  Unfortunately  several  of  his  observa- 
tions have  been  proven  only  for  particular  cases. 

Observing  that  all  nonlinear  systems  of  fluid  dynamics  satisfy 
certain  "conservation  laws",  Lax  states  that  any  hydrodynamic  system  can 
generally  be  brought  into  the  form 

Ut  +  Fx  +  B  =  0,  (4.1) 

where  U  is  a  column  vector  of  unknown  functions,  F  is  a  column  vector 
such  that  F  =  F(x,t,U),  and  B  is  a  vector  coefficient.  U  is  said  to  be 
a  weak  solution  of  eqn.  (4.1)  with  initial  value  $  if  the  integral 
relation 


00  oo 


/       /      W  U  +  WxF  -  WB}  dxdt  +  /     W(x,0)$(x)dx     =     0     (4.2) 

0         -oo  -00 
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holds  for  every  test  vector  W  which  has  continuous  first  derivatives  and 
vanishes  outside  of  some  bounded  region  in  the  x,t-plane.  Equation  (4.2) 
is  obtained  by  multiplying  eqn.  (4.1)  by  W  and  then  integrating  by  parts 
as  follows: 


oo     oo 


Of,,/,/   {WU.  +  WFV  +  WB}  dxdt 
0   -oo    z  x 


OO      00 


/  [WU]  dx ■ -  /   /  UWtdxdt  +  /  [WF]   dt 

-oo       0        0    -oo  0        -oo 


OO      00  00      oo 


/   /  FW  dxdt  +     f      f    WBdxdt  (4.3) 

0     -oo     X  0     -oo 


h  -   '2+I3  +  I5  • 


Now  I-  =  0  since  W  vanishes  outside  some  bounded  region.  Similarly 


^  =  -  /  W(x,0)$(x)dx.  (4.4) 


Hence  we  obtain  eqn.  (4.2). 

Clearly,  the  only  requirement  for  U  to  be  a  weak  solution  is  that 
it  be  continuous  almost  everywhere  so  that  the  (Riemann)  integrals  in 
eqn.  (5.2)  exist.  Hence  weak  solutions  need  not  be  differentiate. 
The  motivation  for  developing  the  notion  of  weak  solutions  is  that  in 
physical  systems  we  are  concerned  with  discontinuous  functions  that 
satisfy  eqn.  (4.1)  almost  everywhere. 

The  Rankine-Hugoniot  equations  now  have  an  interesting  interpretation 
in  terms  of  weak  solutions;  if  U,  and  IL  are  two  genuine  solutions  of 
eqn.  (4.1)  whose  domains  in  the  x,t-plane  are  separated  by  a  smooth  curve, 
the  two  taken  toqether  will  constitute  a  weak  solution  if  and  only  if  the 
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slope  m  of  the  separating  curve  and  the  value  of  U,  and  IL  on  the  curve 
satisfy 

1  (U]  -  U2)  =  F(U1)  -  F(U2)  (4.5) 

Equation  (4.5)  corresponds  precisely  to  eqns.  (2.21)  and  (2.22). 

Havinq  generalized  the  concept  of  a  solution  to  a  differential 
equation,  we  night  expect  some  of  the  properties  of  the  original  concept 
to  be  lost.  For  instance,  initial  values  do  not  in  qeneral  determine  a 
unique  weak  solution.  This  property  raises  immediate  difficulties.  In 
nature  nonviscous  flow  is  described  by  a  unique  weak  solution  to  the 
hvdrodvnamic  equations,  given  an  initial  vector.  Hence  if  our  mathemati- 
cal model  is  a  meaninqful  one,  there  must  be  some  other  Drinciple  that 
uniquely  defines  a  weak  solution.  Lax  offers  some  possibilities,  the 
most  likely  beinq  that  the  weak  solutions  occuring  in  nature  are  limits 
of  viscous  flows.  It  should  be  possible  then  to  determine  some  relation- 
ship between  weak  solutions  and  solutions  to  the  viscous  flow  problem. 
THEOREM: 

Consider  the  nonlinear  parabolic  system 

Ut+Fx+B     -     AUXX  (4.6) 

with  initial  vector  U(x,0)  =  $(x).  Here  AUxx  corresponds  to  a  viscosity 
factor.  Given  that  the  strong  limit  (as  x->0)  of  the  net  of  solution 
functions  U  (x,t)  exists  and  is  equal  to  U(x,t),  then  U(x,t)  is  a  weak 

A 

solution  of  eqn.  (4.1 ) . 
PROOF: 

Consider  an  arbitrary  twice  differentiable  test  vector  W.  Multiplying 
eqn.  (4.6)  by  W  and  integrating  by  parts  yields 


3/ 


00      00 


/    /  (W,UX  +  W  F(UX)  -  WB)  dxdt  +  /  W(x,0)$(x)dx 

Q      -00  -oo 


Xf      f    W  U.  dxdt 

0   -oo    X  A 


(4.7) 


Now  keepina  <£>  and  W  fixed,  and  letting  A-*0,  the  left  side  of  eqn.  (4.7) 
approaches  the  left  side  of  eqn.  (4.2),  and  the  riqht  side  of  eqn.  (4.7) 
approaches  zero.  Hence  U(x,t)  satisfies  eqn.  (4.2)  for  all  twice 
different!' able  test  functions.  Note  that  U  had  to  be  a  strong  limit 
of  U,  (i.e.  //  III,  —  U |  -*-  0  over  any  bounded  region  in  the  x,t-plane) 

A  A 

in  order  that  F(U,)-HF(U).  If  Ih-HJ  only  in  the  weak  sense  (i.e. 
U,(x,t)  U(x,t)  \/  (x,t))there  is  no  quarantee  that  F(U,)  would  converge 
to  F(U).  Now  it  must  be  shown  that  eqn.  (4.2)  is  satisfied  for  any 
arbitrary  test  vector  W  which  is  continuously  once  differentiable.  Lax 
states  that  since  U  satisfies  eqn.  (4.2)  for  all  twice  differentiable 
test  vectors  W,  a  fortiore,  U  satisfies  eqn.  (4.2)  for  all  once  different- 
iable test  vectors.   This  author  disagrees  with  this  argument  since  the 
class  of  once  differentiable  functions  is  laraer  than  the  class  of  twice 
differentiable  functions.  Though  possibly  true,  the  statement  requires 
proof.  Several  approaches  were  attempted  and  proved  to  be  unsuccessful. 
One  apparent  way  to  eliminate  the  difficulty  is  to  restrict  our  attention 
to  just  twice  differentiable  test  vectors  in  our  original  definition  of  a 
weak  solution. 

Having  shown  for  a  particular  case,  that  the  limit  of  viscous  flow 
(as  \-*0)   is  a  weak  solution,  Lax  conjectures  that  all  viscosity  methods 


Lax,  P.  D. ,  "Weak  Solutions  of  Nonlinear  Hyperbolic  Equations  and 
Their  Numerical  Confutation,"  Communications  on  Pure  and  Applied 
Mathematics,  v.  7,  n.  163,  1954. 
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should  converqe  to  the  same  weak  solution.  However  he  proposes  a 
different  limitinq  process,  namely  a  special  finite  difference  scheme, 
which  is  independent  of  the  viscosity  of  the  fluid.  This  aoDroach  has 
a  significant  advantaqe  over  the  viscosity  aoproach  in  that  it  accurately 
represents  the  actual  behavior  of  the  fluid  in  e\/ery   region  of  the  flow. 
Lax's  method  replaces  the  differentiations  by  finite-difference  operations 
accordino  to  the  scheme 

fx(x,t)+  ^x"  Cf(x+Ax,t)  -  f(x-Ax,t)] 

(4.8) 

ft(x,t)-  ^  [f(x,t+At)  -  nx+AX,t)  -  f(x-AX,t)] 

Finally,  substantial  numerical  evidence  supports  Lax's  conjecture 
that  when  the  above  finite  difference  scheme  is  applied  to  any  single 
homogeneous  first  order  conservation  law 

ut  +  [f(u)]x  =  0     f "  <  0  (4.9) 

with  U(x,0)  =  <f>(x),  the  solution  U(x,t,Ax)  approaches  the  same  weak 
solution  generated  by  the  viscosity  methods.  The  solution  is  given  by 

u(x,t)  =  g  (-^-)  (4.10) 

where  y   =  y  (x,t)  maximizes 
Jo  Jo 

/  <o(s)ds  +  tG(^)  (4.11) 

0  z 


and 


f'[g(s)]  =  s  ;   G'(s)  =  g  (4.12) 


In  most  physical  situations  U(x,t)  is  a  piecewise  differentiable  function, 
Then  it  is  an  easy  matter  to  show  U(x,t)  is  a  weak  solution  to  eqn.  (4.9) 
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since  verification  of  eqn.  (4.2)  may  be  avoided.  All  that  must  be  shown 
is  that  U(x,t)  satisfies  eqn.  (4.9),  wherever  it  has  well-defined  first 
derivatives.  Since  the  first  derivatives  of  U(x,t)  are  undefined  at  most 
on  a  set  of  measure  zero,  and  U(x,t)  satisfies  eqn.  (4.9)  everywhere 
else,  U(x,t)  must  by  necessity  satisfy  the  intearal  relation  (4.2) 
everywhere.  Hence  (4.10)  is  a  weak  solution  of  eqn.  (4.9).  The  intearal 
eqn.  (4.11)  is  necessary  to  uniquely  define  the  weak  solution. 
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V.  CONCLUSIONS 

The  original  intention  of  this  paper  was  to  analyze  difficulties  in 
a  numerical  system  called  "PUFF"  [Ref.  2],  "PUFF"  is  an  attempt  to 
numerically  compute  the  reaction  of  a  multi-layered  medium  to  violent 
shocks.  The  Puff  code  employs  Von  Neumann  and  Richtmyer's  artificial 
viscosity  term.  It  has  been  shown  that  there  were  very  large  discrep- 
ancies between  the  "PUFF"  solution  and  the  classical  solution  inside 
the  shock  regions  [Ref.  1].  The  reason  for  these  errors  should  now  be 
apparent.  Von  Neumann  and  Richtmyer's  artificial  viscosity  term  had 
to  comply  with  certain  requirements.  However  these  requirements  were 
to  a  large  degree  independent  of  actual  physical  considerations  and, 
hence,  were  not  sufficiently  restrictive.  Consequently,  in  shock 
regions  where  viscosity  is  significant,  the  artiftcfal  term  may  not 
accurately  describe  the  actual  fluid  flow.  This  difficulty  becomes 
critical  in  a  multi-layered  medium  since  the  number  of  shock  waves  is 
increased  when  the  original  shock  reflects  from  several  boundaries. 
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